function [e_Y_cpi,m_e_mj_cpi,cov_e_mj_cpi,check_mj_cpi,...
    e_Y_doubledef,m_e_mj_doubledef,cov_e_mj_doubledef,check_mj_doubledef,...
    e_Y_def,m_e_mj_def,cov_e_mj_def,check_mj_def]...
    = elasticity_decomp_sectorFRCES_fun(delta,france,s1,s2,VA_hat_mj,P_hat_n,...
    RGDP_cpi_hat_n,VA_mj_0,GDP_deflator_n,RGDP_def_hat_n,DoubleDeflation_deflator_n,RGDP_doubledef_hat_n)


global shockT alphaT rhoT lambdaT etaT rho sigma lambda eta varphi ...
    tol maxit D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    N J FI_J france countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO sector_algorithm ... 
    sL_n sPi_n sD_n alpha_WIOD_j alpha_WIOD_francej labour_share_PWT

%% Calculation using CPI

% Calculates the elasticity decomposition at the sector (ij) level for
% France only, so will be based on data aggregated from firm-level data for
% the firm-level model

% Decomposition of aggregate real GDP elasticity into sector mean and
% covariance components

% Note: N = 15 for France

VV_mj = VA_hat_mj(s1:s2,:);
real_VA_hat_mj = VV_mj/P_hat_n(france);

e_mj_cpi = (real_VA_hat_mj-1)/delta;
e_Y_cpi = (RGDP_cpi_hat_n(france)-1)/delta;

m_e_mj_cpi = nanmean(e_mj_cpi,1);

VV_mj = VA_mj_0(s1:s2,:);
w_mj = VV_mj/sum(VV_mj);

w_mj_bar = size(w_mj,1)-1;  % N - 1 correction for small small as in COV
clear  VV_mj;

C = nancov(w_mj*w_mj_bar,e_mj_cpi);
cov_e_mj_cpi = C(1,2);

check_mj_cpi = e_Y_cpi-m_e_mj_cpi-cov_e_mj_cpi;

%% Calculation using deflator

VV_mj = VA_hat_mj(s1:s2,:);
real_VA_hat_mj_def = VV_mj/GDP_deflator_n(france);

e_mj_def = (real_VA_hat_mj_def-1)/delta;
e_Y_def = (RGDP_def_hat_n(france)-1)/delta;

m_e_mj_def = nanmean(e_mj_def,1);

VV_mj = VA_mj_0(s1:s2,:);
w_mj = VV_mj/sum(VV_mj);

w_mj_bar = size(w_mj,1)-1;  % N - 1 correction for small small as in COV
clear  VV_mj;

C_deflator = nancov(w_mj*w_mj_bar,e_mj_def);
cov_e_mj_def = C_deflator(1,2);

check_mj_def = e_Y_def-m_e_mj_def-cov_e_mj_def;

%% Calculation using double deflator

VV_mj = VA_hat_mj(s1:s2,:);
real_VA_hat_mj_doubledef = VV_mj/DoubleDeflation_deflator_n(france);

e_mj_doubledef = (real_VA_hat_mj_doubledef-1)/delta;
e_Y_doubledef = (RGDP_doubledef_hat_n(france)-1)/delta;

m_e_mj_doubledef = nanmean(e_mj_doubledef,1);

VV_mj = VA_mj_0(s1:s2,:);
w_mj = VV_mj/sum(VV_mj);

w_mj_bar = size(w_mj,1)-1;  % N - 1 correction for small small as in COV
clear  VV_mj;

C_doubledeflator = nancov(w_mj*w_mj_bar,e_mj_doubledef);
cov_e_mj_doubledef = C_doubledeflator(1,2);

check_mj_doubledef = e_Y_doubledef-m_e_mj_doubledef-cov_e_mj_doubledef;

end

